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^ ; ABSTRACT 

We calculate the rate of hydrogen burning for neutron stars (NSs) with hydrogen atmospheres and 
an underlying reservoir of nuclei capable of proton capture. This burning occurs in the exponentially 
suppressed diffusive tail of H that extends to the hotter depths of the envelope where protons are rapidly 
captured. This process, which we call diffusive nuclear burning (DNB), can change the H abundance at 
the NS photosphere on timescales as short as 10 2-4 years. In the absence of diffusion, the hydrogen at 
the photosphere (where T m 10 6 K and p ~ 0.1gcm~ 3 ) would last for far longer than a Hubble time. 
Our work impacts the understanding of the evolution of surface abundances of isolated NSs, which is 
important to their thermal spectrum and their effective temperature-core temperature relation. In this 
paper, we calculate the rate of H burning when the overall consumption rate is controlled by the nuclear 
timescales, rather than diffusion timescales. The immediate application is for H burning on millisecond 
radio pulsars and in quiescence for the accreting NS Cen X-4. We will apply this work to young radio 
pulsars and magnetars once we have incorporated the effects of strong B > 10 12 G magnetic fields. 

Subject headings: conduction - diffusion - stars: abundances, interiors - stars: neutron - 
' X-rays:binaries 
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1. INTRODUCTION 



The increasing number of observations of isolated neutron star (NS) atmospheres (see Pavlov, Zavlin & Sanwal 2002 
I ' for a recent review) and resulting constraints on their surface composition has highlighted the need to consider the role of 
Q , nuclear physics during their cooling phase. The initial composition of the outer layer is neither known nor constrained by 
(— > ■ the theory of supernova explosions since the amount of matter needed to affect the outer envelope is small (~ 10~ 20 M Q 
^ | for the photosphere). This much material can easily fallback and contaminate the outer atmosphere (Woosley & Weaver 
1995). Even if the fallback material consists of heavy elements, spallation (see Bildsten, Salpeter & Wasserman 1992) will 
occur at some stage in the fallback, creating a plethora of lighter elements that rapidly gravitationally separate. Current 
evidence suggests that young pulsars (< 10 4 ~ 5 yrs) possess photospheres of either hydrogen or helium, whereas older 
pulsars (> 10 5 yrs) appear to have envelopes of heavier, more uniformly opaque elements (Yakovlev, Kaminker & Gnedin 
Oj 2001). This may indicate an evolution of the outer envelope on the timescale of 10 4_5 yrs, something we would like to 
understand better. 

As part of our ongoing work on this problem, we present here a mechanism of nuclear processing of the outer envelope 
which we call diffusive nuclear burning (DNB). The basic picture is simple and was alluded to by Chiu & Salpeter (1964) 
as a mode of nuclear burning. Rosen (1968) carried out the initial calculations of the effect for surface temperatures in 
the range of (5 — 10) x 10 6 K. Their work was in the limit where the diffusion time to larger depths was the rate-limiting 
step. 

In this paper, we will consider much lower values for the surface temperature, where diffusion is not the limiting step. 
Consider a NS envelope that consists of hydrogen above a layer of carbon (or any other proton capturing material) in 
diffusive equilibrium (see illustrative Figure 1). For simplicity, we do not calculate the effect of an intervening helium 
layer in this paper, though we estimate its possible impact in the conclusions. As long as there has been adequate time 
to reach a diffusive equilibrium, the separation between the hydrogen and carbon is not strict and a diffusive tail of 
hydrogen penetrates deep into the carbon layer. The temperature rises with depth and can increase by 1-2 orders of 
magnitude 10 meters beneath the photosphere. At this location, the lifetime of a proton is very short. The convolution of 
the increasing nuclear burning rate and the decreasing hydrogen abundance creates a burning layer at this depth, where 
hydrogen burning peaks. As hydrogen burns, the diffusive tail is driven slightly away from diffusive equilibrium and a 
diffusive current is set up, «7dnb- We will always work in the limit where Jdnb is far less than the maximum current 
that can be supported due to diffusive processes so that diffusive equilibrium remains an excellent approximation. 

We find that, even when the photospheric temperatures and densities are so low that the local nuclear lifetime of a 
proton is in excess of a Hubble time, DNB can proceed on a much faster timescale. In this paper we make some simplifying 
assumptions that prohibits a direct application to radio pulsars, and thus the extension of the physics to radio pulsars and 
magnetars will be left to future work. For example, we require a magnetic field low enough so that there is a negligible 
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effect on the temperature profile; B < 10 9 G (Ventura & Potekhin 2002; Potekhin & Yakovlev 2001). However, our work 
is applicable to some X-ray transients in quiescence, where the outer envelope consists of pure hydrogen, even though the 
depth of the hydrogen layer is unknown. 

This paper is structured as follows. In § 2 we present the basic equations and microphysics governing the thermal and 
compositional structure of a NS envelope in diffusive equilibrium. The derivation of the electric field strength from the 
envelope structure equations is in § 3. Using the derived electric field strength, we numerically solve the thermal and 
compositional structure of the NS envelope. The analytic solution to the NS envelope structure is given in § 4. In § 5 
we discuss the relevant hydrogen burning reactions that drive the NS envelope evolution, present numerical and analytic 
solutions for Jdnb, and apply our theory to the X-ray transients, specifically Cen X-4. We summarize our findings in § 6. 

2. ATMOSPHERIC STRUCTURE AND MICROPHYSICS 

We only consider the top 10 4 cm of the envelope, where the density is < 10 10 g cm~ 3 . Since the thickness of the envelope 
is much smaller than the NS radius (R w 10 km), we presume a plane parallel atmosphere with a constant downward 
gravitational acceleration, g = GM/R 2 . We neglect all relativistic corrections. The structure of the NS envelope is 
determined by the equations of hydrostatic balance, 

dPi 

= —Hi (Airripg - ZieE) , (1) 

dP e 

= -n e (m e g + eE) , (2) 

where Pi, rij, Ai, Zi are the pressure, number density, atomic number and charge of the z'th ion species and E is the 
upward pointing electric field. The thermal structure is determined by the heat diffusion equation for a constant flux, 
acT^/4, where T e is the effective temperature, 

dT_ 3k P 4 

dr 16T3 e ' w 
where k is the opacity. Finally, we demand charge neutrality, 

n e = ^2 l Z i n i . (4) 

i 

Summing the pressure equations with the charge neutrality constraint recovers the familiar form of hydrostatic balance, 
dP/dr = —pg. 

In the outermost layers of a NS envelope, the nuclei are fully ionized and the radiative opacities are set by free-free 
absorption and Thomson scattering, 

n e a Th 0.4 cm 2 

K T h = ~ , (5) 

P Me g 

where p e = A/Z is the electron mean molecular weight. We calculate the free-free opacity by adopting the formalism 
used in Schatz et al. (1999), 



k„ =7.53xl0 6 ^^^^.g ff (Z 4 ,T,n e )^, (6) 



where T§ = T/10 6 K, p$ = /9/10 5 gcm~ 3 , Xi is the mass fraction, and gs is the Gaunt factor. For our numerical 
calculations, we calculate the Gaunt factor using the Schatz ct al. (1999) fitting formula [given by their equation (A2)]. 
For our analytic calculations we set gs = 1. 

At very high densities, the dominant mode of heat transport is conduction by degenerate electrons. For strongly 
coupled Coulomb plasmas (SCCP), the electron thermal conductivity is determined by electron-ion scattering as given by 
the Wiedemann- Franz law (Ventura & Potekhin 2002; Brown, Bildsten & Chang 2002). 

where m* = cf/c 2 is the effective electron mass, €f is the electron Fermi energy, and tk is the effective relaxation time. 
For weak magnetic fields (B < 10 9 G), we use the thermal conductivity of electrons in SCCP given by Baiko et al. (1998) 
with analytic formulae given by Potekhin ct al. (1999). In the liquid metal approximation, Brown et al. (2002) showed 
that this formalism agrees with the calculations of Urpin & Yakovlev (1980) and Itoh et al. (1983). 

The transition of the envelope from primarily radiative to primarily conductive takes place over a narrow region called 
the sensitivity strip (Ventura & Potekhin 2002). Heat transport by radiation and conduction is of the same order in this 
region. The sensitivity strip also plays the dominant role in determining the temperature profile of the NS envelope (as 
discussed in § 4, also see Ventura & Potekhin 2002). In the sensitivity strip, the SCCP is always in the liquid metal 
regime, whose conductivity is fairly well modeled (Brown et al. 2002; Ventura & Potekhin 2002). Thus, we are confident 
in our results for the thermal profile. 
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The electron equation of state transitions from an ideal gas EOS to a degenerate EOS. To model this transition in 
our numerical solutions, we have used the Paczynski (1983) equation of state. In our analytic calculations, we solve our 
envelope for the two different limits and introduce a free parameter ^ to join them at a point defined by 

, (8) 

e,nondcg 

where -P c ,dc g and P c , n0 ndcg are the nonrelativistic degenerate and nondegenerate electron equations of state respectively. 
For this calculation we take ^ = 1.5, which gives a reasonable approximation to the full numerical solution in the region 
of the envelope where DNB occurs. 

3. ELEMENTAL DISTRIBUTION AND ELECTRIC FIELD STRENGTH 

The composition of the envelope in diffusive equilibrium is found by simultaneously solving the equations of hydrostatic 
balance for each species and the flux equation. This approximation is excellent once enough time has passed for the 
diffusive equilibrium to be established (Brown et al. 2002). 

We begin by calculating the local electric field. In early studies of white dwarfs and NSs (Rosen 1968; Fontaine & 
Michaud 1979), the electric field was presumed to be that for an ideal gas of ions and electrons, 

eE = -^—^m p g. (9) 

However, in the degenerate case where P e ~ P, the field is 

eE = ~z m p g - ( 10 ) 

These results are only valid in certain limits, where the plasma consists of one dominant species and where the gas is either 
very degenerate or very non-degenerate. For multi-ion plasmas and for partial degeneracy, these simple approximations 
break down. For a more general solution, we calculate the electric field from the ion and electron equations of hydrostatic 
balance (cq. [1]), and the flux equation (eq. [3]). We presume an ideal gas equation of state for the ions and an arbitary 
electron equation of state. Finally we drop the m e g term. 

Let us first consider the isothermal case. We first expand dP e /dr, 

dn e _ / dP e \ 
dr \ dn e 

and then use the ideal gas equation of state for the ions: 



n e eE, (11) 



^ = -^ {Ai m p9 -eZ i E). (12) 

Multiplying each ion equation (eq. [12]) by Zi and subtracting their sum from the electron equation (eq. [11]) and 
presuming charge neutrality (eq. [4]) gives 

eE= ™ P 9En t A tZl 

E^Zf +n e k B T(dP e /dn e ) 

For the nonisothermal case, we get an additional term in the expansion of dP e /dr and dPi/dr, e.g. 

and the electric field becomes 

cE _ J2 nA (A imp g + k B (dT/dr)) - (dPJdT) (dPjdn.y 1 (dT/dr) 

Y,n % Z^ + n e k B T{dP e /dn e )- 1 

This result agrees with the results of Macdonald, Hernanz & Jose (1998) and Althaus & Benvenuto (2000), but we have 
decoupled our solution of the electric field from the drift velocity and diffusion coefficient. Hence, the electric field is 
simply a function of local parameters. For the purposes of numerical calculation we use the Paczynski (1983) electron 
equation of state. For our analytic calculations, we use either the ideal gas or the nonrelativistic degenerate electron 
equation of state depending on the local conditions of the plasma. 

At this point, it is helpful to illustrate a solution for the particular case of a H/ 12 C envelope. In Figure 2, we plot the 
electric field of the layer as a function of column depth, y = P/g. We integrate from the photosphere with a background 
of H. At a certain depth we introduce a low carbon abundance and continue integrating inward. As expected, the electric 
field approximates that of the value given by equation (9) at the extreme limits where one ion dominates over the other. 
The slight discrepancy between the value give by equation (9) at the extreme limits is due to degeneracy effects from the 
Paczynski (1983) fit of the electron equation of state. The flux and structure equations can now be solved self consistently 
to model any arbitrary envelope in diffusive equilibrium. In Figure 3, we plot the solution for a H/ 12 C envelope. We plot 
the composition profile in terms of the concentration of hydrogen, f B = nn/n-t, where n t is the total number of ions. 
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4. APPROXIMATE SOLUTIONS FOR TEMPERATURE AND COMPOSITION 

The temperature and composition profiles can be represented well with approximate analytic solutions that illuminate 
the physics. Hernquist & Applegate (1984) and Ventura & Potekhin (2002) have presented analytic thermal calculations of 
the NS envelope, but their solutions are too complex for our analytic estimates. Hence, we present simplified approximate 
solutions in § 4.1. For an envelope in diffusive equilibrium, Fontaine & Michaud (1979) have presented an excellent 
approximate solution for the composition profile. In § 4.2, we extend their technique to the degenerate regime. 

4.1. Analytic Temperature Profile 

Roughly speaking, the NS envelope can be divided into three zones, a radiative outer zone, a sensitivity strip (where 
opacity changes from primarily radiative to conductive) , and a conductive isothermal interior (also see Ventura & Potekhin 
2002). In the radiative, nondegenerate outer zone, the opacity is mostly determined by free- free absorption. With an ideal 
gas equation of state, yg — pksT/ fim p , where jj, = A/ (Z + 1) is the mean molecular weight of the plasma, the solution to 
the heat diffusion equation and hydrostatic balance is 

T(y) = 1.58 x 10 V/ 17 (fiZg 14 T^~ 2 ) 2/17 K, (16) 
where for brevity we have written gi4 = g/(l0 14 cms -2 ) and T e $ = X" e /(10 6 K). For a radiative nondegenerate envelope 
consisting of multiple layers of ions, the solution is modified via the introduction of constants of integration at the 
boundaries. If the boundaries are sufficiently far away in pressure from the sensitivity strip, the temperature profile of 
the underlying layers will approximate that of a pure component. 

In the sensitivity strip, which is defined by K r ad(j/ss) = ^cond(j/ss), the opacity transitions from primarily radiative to 
primarily conductive. Here we follow the formalism of Ventura & Potekhin (2002). To find the sensitivity strip, we take 
the approximate expression for conductivity (Ventura & Potekhin 2002), 

g c = 2.3xl0 15 ^ X \ 6rg ° (17) 
AZ 1 + X 2 cmsK' v ; 

where \ — PF/ m e c ~ (Pe/pe) 1 ^ 3 is the relativity parameter, [i e = A/Z, and A is the Coulomb logarithm. We set A = 1 

for simplicity's sake and 1 + \ 2 ~ 1 since the relativity parameter is small in the nonrelativistic regime. Equating the two 

opacities, equations (6) and (17), we find that the sensitivity strip is defined by the condition, 

A l/3 T 17/6 2/3 

V = 11-2 r,g Me -V (18) 
Z * gull enr 

Taking this condition and inserting the solution for the radiative zone, we obtain expressions for the column depth of the 
sensitivity strip, y ss , 

and temperature, 



y ss = 6.86 x 10 4 ^% (19) 



/ v3rpl2 \ 2 / 17 

T ss = 2.2 x 10 7 -g-f K, (20) 

as a function of effective temperature, gravity, and composition. Beyond the sensitivity strip (y > y ss ), we solve the 
constant flux equation using the conductive opacity, 

T— — 2.2 x 10 16 ^^f^-, (21) 
ay p l 

when 1 + x 2 ~ 1- Integration of this equation is relatively straightforward given an equation of state. There are two cases 
that we need to consider. For the nondegenerate regime, the solution for y > y ss is 



T = T ss exp 



-1.62 x 10 4 ^|^ Or 1 -y?) 
9uP 



K, (22) 



where we have joined it with the radiative solution at the sensitivity strip. 

The electrons become more degenerate as the column increases. We define the degenerate boundary with our analytic 
degeneracy condition (eq. [8]), which gives the condition 

T 5/2 5/2 5/3 

y = ^ 6 5/ f (23) 
gufi a/ cm z 

where we have assumed that ydog 3> Vss- Using our solution for T beyond the sensitivity strip (eq. [22]), we find a system 
of equations for Td eg and j/dcg, which are solved numerically. 
The solution beyond the degenerate boundary is 



Ti„ = 1.37 x 10 16 T, 4 



Z 



(»i /6 -»- 1/B )K a . (24) 



L deg ± -°' * ± e& 6/5 

9lA Pe 

Taking the limit of equation (24) as y — > oo gives an analytic expression for the core temperature, T c , as a function of 
T e , which reproduces previous scalings (Gudmundsson, Pethick & Epstein 1983; Potekhin, Chabrier & Yakovlev 1997). 
Ventura & Potekhin (2002) present a more generalized solution in terms of the relativity parameter \- 

In Figure 3 we present the approximate and numerical solutions for the temperature and composition which we discuss 
in the next section. For the thermal profile, the agreement is excellent (to within 10%). It also agrees well with the 
analytic solution given by Ventura & Potekhin (2002). 
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4.2. Analytic Composition Profiles 

We find analytic solutions for the composition when one ion species is dominant and the electrons are either nondegen- 
erate or very degenerate. For the nondegenerate case, we first derive the result of Fontaine & Michaud (1979). We start 
with the hydrostatic balance equations for the ions, where "1" is the background and "2" is the trace, 

dPi 
dr 
dP 2 
dr 

For the ions, Pi = mk B T. Dividing each pressure equation above by their respective Pi and subtracting "1" from "2", we 
have: 

1 d\n(P 2 /P 1 ) 



= -m (A±m p g - Z-i.eE) , 
= -n 2 (A 2 m p g - Z 2 eE) . 



(25) 
(26) 



k B T- 

which in the trace limit (n 2 <C n\) gives 



dr 



{A 2 -A 1 )m p g+{Z 2 -Z 1 )eE, 



hi'* 
Pi 



ln(^ 

m 



In 



n 2 



rii + n 2 



In/: 



2- 



(27) 



(28) 



where f 2 is number fraction of the trace. 

In the nondegenerate regime, the electric field is given by equation (9). Using this result and using the trace approxi- 
mation so that A « Ai and Z ss Z\, equation (27) becomes 



Using hydrostatic balance dP/dr 



d\nf 2 

dr 
-pg « 



A\m p g 
(Zi + l)k B T 
-Aim p gni and P 



M 
Ai 



(Zi + 1) - Z 2 - 1 



(29) 



din / 2 



^2 



(^i 



dlnP 

The concentration of a trace species is a power law of pressure with an exponent, 8 
applicable because it gives the concentration in terms of a pressure contrast, 

MP) 



[Z\ + l)n\kT, we have the elegant result, 
l)-Z 2 -l . (30) 
A 2 (Zi + l)/Ai-Z 2 -l. It is widely 



Pbj 



(31) 



where f 2 (Pb) is the concentration at an arbitrary boundary P^. 

We now extend this result to the highly degenerate regime, where the electron pressure dominates over ion pressure, 
and the electric field is 



eE = 



m p g 

He 



A 



-m p g. 



Placing this into equation (27) gives 



din / 2 

dP p(P,T)k B T(P) 



AiTUp 



A x 



?1 
Zi 



(32) 



(33) 



We solve for p{P, T) via the electron equation of state and T(P) from the flux equation. This solution is best when the 
electrons are highly degenerate. However, this is also the regime where the temperature begins to be isothermal (Ventura 
& Potekhin 2002), allowing us to replace T(P) by Td eg , where Td eg is the temperature at the degenerate boundary, giving 



A\m p 



rfln/a = 

dP Pd e S (P/Pde e ) lh k B T d 



Ai 



Zn 
Zi 



(34) 



where 7 is the polytropic index of the degenerate electron equation of state, Pap 7 . The condition of degeneracy (eq. 
[8]) determines Pdeg and Pdcg- 
The electrons are nonrelativistic at the degenerate boundary so 7 = 5/3 and we integrate (eq. [34]) to get 



H = /2,deg,0 exp <^ -T) 



PdogfcijPdog 



Z2 

Zi 



Ai 



p 

~p&~ ci 



- 1 



(35) 



where /2,dc g ,o is the concentration where the envelope becomes degenerate and 77 = 1 — - = 0.4. Note that 77 is positive 
for all degenerate equations of state: the concentration always decreases exponentially with increasing pressure. Our 
solution reproduces the Boltzmann solution for a particle experiencing a upward pointing force F in the isothermal limit, 
f 2 cc exp(FrfkT) since r oc P/p oc pi. Thus, in the degenerate regime the concentration falls faster than the power 
law relation in the nondegenerate regime. In Figure 3, we compare the numerical solution with our analytic calculation. 
The agreement is poor for higher densities due to the increasing importance of relativistic equation of state, but gives 



reasonable agreement in the burning layer ybu 



10 6 
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5. HYDROGEN BURNING IN DIFFUSIVE EQUILIBRIUM 

The relevant nuclear processes that consume hydrogen depend on the temperature and composition of the underlying 
matter. We consider the case of a NS which has a pure hydrogen layer on top of a layer of proton capturing material (i.e. 
carbon, nitrogen, oxygen, etc) in diffusive equilibrium. We have not yet calculated the impact of an intervening He layer, 
but we think it will be small as long as it does not penetrate to a depth greater than the burning layer. 

For simplicity, we first consider carbon. For hydrogen on carbon, there are only two possible nuclear reactions: 
12 C(p,7) 13 N and p(p, e + + v e )T). The D produced in the p-p capture undergoes D(p, 7) 3 He which in general is the 
cndpoint of this process. Since we always have a region where T > 10 7 K, the reaction rate for p-p capture is very slow 
compared to the proton capture rate onto carbon. The 13 N decays to 13 C which has A/Z = 13/6 « 2.167, greater than 
the local A/Z = 2 for pure carbon. The 13 C will thus sink through the outer layers of the NS, and will not reside long 
enough in the burning layer to facilitate a catalytic cycle. Hence, the complete exhaustion of H by this process requires 
an excess of 12 C compared to H. 

5.1. Diffusive Nuclear Burning 

We calculate the burning rate presuming that the burning time is slow compared to the time for hydrogen to diffuse 
down through the carbon to the burning layer. We express this hierarchy of diffusion time to nuclear burning time by 
studying the H equation of continuity, 

dn H n H 

-— — + V-J H = , 36 

dt t h 

where t h = (<jv)~ 1 n~ c 1 is the local lifetime of H to 12 C capture. The condition of diffusive equilibrium is drin/dt = 
or more specifically that dun / dt changes only on the timescale associated with depletion of the hydrogen column, which 
is long compared to both the local nuclear burning time tr and the diffusion time and therefore is dropped. With 
the condition of steady state J h and Jdnb are equivalent. Writing the result in one dimension for J h = Jdnb = 

Vdnn/dz — n#t>d r , we have 

d 2 n H dn H n H 
oz z oz t h 

, where V is the diffusion coefficient and Wdr is the drift velocity. We presume that V and Wdr change slowly with z. Our 
presumption that nuclear burning does not affect diffusive equilibrium implies that locally the timescale associated with 
nuclear burning is much longer than the timescale associated with diffusion to the burning layer, or th 3> Tion = I 2 I'D, 
where I is the ion scale height for Hin the 12 C layer. Therefore, taking T> from Brown et al. (2002), we have 

A^gl^Ai-Zif 

where Ai and Zi are the atomic number and charge of the background proton capturing element. 
For proton capturing reactions, the proton lifetime (Clayton 1983) is 



r ion = 800 , n A /, ' 2 s, (38) 



I =^-° ,a (ii^)¥^ /a -(-^- ,/a )-- 



(39) 



where Xi and Ai is the mass fraction and atomic number of the proton capturing element and So and B are parameters 
determined by the proton capturing element. We derive an analytic expression for the condition of diffusive equilibrium 
by expanding equation (39) about T§ = 40, 



-14 



t h = 1.32 x 10 5 P 5 1 (^J s. (40) 

Comparing equations (38) and (40), our condition for carbon is 

T 6 < 45.5p 5 - 011 , (41) 

for #14 = 2. 

If the conditions for DNB are fulfilled for a sufficiently large pressure range, the total burning rate of hydrogen, defined 
as (h = Jdnb (H/ 12 C boundary) is 

a Vh f n Hm p 

Ch = = / — -, * L -^n d *i 42 ) 

Tcoi J TH(nH,nc,T) 

where yn is the integrated column of hydrogen, r co i is the characteristic time for that column to be consumed and nn 
and nc are the local number density of hydrogen and carbon respectively. This equation is incorporated into our flux 
and structure equations and solved simultaneously, yielding t co \ as a function of our stellar properties and yn- We plot 
dC,H I 'd\g w y for an envelope with yjj = 100gcm~ 2 in Figure 4. In the bottom graph we show the local concentration 
and temperature. Due to the inverse dependence of the local hydrogen burning rate on temperature and concentration, 
the burning layer is concentrated in a narrow pressure range at a depth where the H abundance is small. The diffusion 
timescale to this depth is T lon = 10 4 s, which is much shorter than the nuclear timescale of 1.5 x 10 7 s. This observation is 
the essence of DNB, namely, that practically all the burning occurs in the exponentially suppressed diffusive hydrogen tail. 
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In Figure 5, we plot the characteristic time, r co i, and characteristic mass burning rate, Mdnb = ^R 2 Dh /t co \, for 
<7i4 = 2 and fixed core temperature, T c . For a central temperature of 5 x 10 7 K, our calculations are valid only up to 
3 x 10 7 gcm~ 2 . Above this column, the assumption of constant flux breaks down since the energy released from DNB is 
comparable to the flux. 

Two results are obvious from this figure. First the relation between lifetime (and hence mass burning rate) and column 
is related via a simple power law between t co \ and yn- We derive this power law in § 5.2. Secondly, the lifetime of the 
entire envelope is set by the lifetime of the photosphere. 

The location of the burning layer must satisfy our derived diffusive equilibrium condition, (cq. [41]) for our calculation 
to be valid. In Figure 6, we plot various models against our condition for DNB in diffusive equilibrium (eq. [41]). The 
burning layers is the region between the two vertical lines in each model. For mo dels with T e > 10 6 K (or T c > 6 x 10 7 K), 
our present calculation is not valid as DNB will occur in the slow diffusion limit. We will expand on this point and study 
the slow diffusion limit in a future paper. 

5.2. Approximate Solution 

The competition between the falling H abundance and rising temperature with depth forces much of the hydrogen 
burning to occur in a narrow zone. As a result, we can solve the burning integral (eq. [42]) over the burning layer by 
approximating it with the method of steepest descents. 

Since DNB occurs deep below the hydrogen layer, Jj« 1, equation (42) then becomes 

<„ = M _ 2.45 x io-4- (tb&s;) / wW'exp (-but'") iVt)jn r' (43) 

where Jh is the local number concentration of hydrogen, which decreases with increasing y. Note that we have made a 
change of variables of dz = yp~ 1 dln(y). 

This integral should be evaluated in both the nondegenerate, radiative regime and the degenerate, conductive regime. 
However, since the majority of the burning occurs in the degenerate, conductive regime, we will evaluate the integral there 
by the method of steepest descents (Clayton 1983), 

/ 2~7T 

g(x) exp(-h(x))dx « g(x Q ) exp(-/i(x )M h „^y ( 44 ) 

where g{x) is a slowly varying function function of x, h(x) is a peaked function, and xo is the solution of h'(xo) = 0, the 
location of the peak of the burning rate. In order to carry out this approximation for equation (43), we make the following 
identifications, 

exp(-/i(x)) -> ypf H cxp (-BT^ 1/3 ^ , 

g(x)^T- 2/3 . (45) 



/ 



L 6 

For the degenerate, conductive regime, the peaked function is 



-h My)} = 8 - ln(y) BT^ 3 - ^ m f^ ( A. _ _L 



y 

J/dej 



■> ' IKh -J-'L- '-de- 

The last term comes from the analytic solution for the composition profile (eq. [35]), where the background is the proton 
capturing clement and the trace is hydrogen. The temperature, T, in the degenerate regime is given by equation (24). 
Since the burning peak yo is defined by a transcendental equation, /i'(ln(t/o)) = 0, we solve it numerically and then 
calculate h"(ln(yo)). For a total column of hydrogen of y# = 100gcm~ 2 , T e = 8 x 10 5 K, and 514 = 2, we find by the 
method of steepest descents, t co \ = 468 yrs, which compares well with the numerical answer of t co \ — 428 yrs. The burning 
layer is centered around yo = 7.2 x 10 5 gcm~ 2 with a local temperature of T(y — y ) = 2.6 x 10 7 K and central core 
temperature T c = 3.3 x 10 7 K. 

Power law scalings of r co i with yu , T e and g can also be determined from this integral. The simplest power law scaling 
involves yn, which we determine as follows. The rate, j/h/t co i, is directly proportional to fjj as in equation (43). From 
equation (31) and (35), Jh a /H,do g ,o = (ydog/ys) 5 , where /H,dc g ,o is the hydrogen number fraction at the degenerate 
boundary and 5 = A^{Z\ + 1)/A\ — Z2 — 1 = —17/12 for H on 12 C. Therefore, t co \ oc y]^ S , which is confirmed from 
comparisons to the numerical results. This scaling is accurate as long as yn < J/dcg- 

The scaling of r co i against T e6 and g\4 can be derived from a power-law expansion of burning integral (eq. [43] ) in terms 
of T. The expansion of the burning integral (eq. [43]) gives 

— « fHTo 2 6 /3 cxp (-BT j/ 3 ) , (47) 

where T 0i 6 = T /10 6 and T is the temperature of the burning layer. Since we are in the degenerate conductive regime, 
T ~T C . The location of the burning region, y and po, scale weakly with T compared to the temperature of the burning 
region, T . Taking the analytic solution to the temperature profile (eq. [24]) and degeneracy point (eq. [23]), the power 

law dependence between T c , T e6 and g u is T oc T c cc T , e 2 6 .g 14 1 ^ 2 . Therefore, the scaling of fu in 12 C goes as, 



f H oc exp 



Pdcgfcs?0 



( P 



- 1 



(48) 
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Expanding this equation at T = 10 8 K gives f H oc (T /10 8 ) - 475 g°4 285 oc T e ° 6 95 5 ° 4 0475 . Expanding the exponential in 
equation (47) for carbon (B = 136.96) gives exp (-£T~ 1/3 ) oc T c 9 - 8 oc T^- 7 g^ 9 . Thus, the scaling for the rate is 

Vh/tcoI oc Tgg ' 3 5i4 4 ' 52 . Hence, the dominant scaling comes from the exponential dependence on temperature of the 
nuclear burning rate with small corrections from the other factors. Putting all the scalings together with our analytic 
calculation, we find 

In Figure 7, we compare the scaling found numerically with this analytic result and find reasonable agreement. The 
agreements for the scalings of yn and gu and the numerical results are good. However, the agreement for the scaling of 
T e and numerical results diverge for low T e . This is due to the changing power law with respect to expanding around 
different values of T c . Expanding around a value of T c appropriate for these lower values of T e resolves this discrepancy. 

These scalings can be calculated for other proton capturing elements. The scaling of t co \ and yn remains the same 
with the exponent 1 — where 5i takes on different values for the differing backgrounds (Zj, Ai). We can use the same 
strategy as earlier, expanding our results in Tq and inserting the scalings for T e g and 314 afterwards. Hence, fjj has the 
same scaling as before, fjj cc (To/10 8 ) a475 g°^ 285 . However, the dominant scaling from the exponential dependence of 
temperature on the burning rate changes. Hence we choose to ignore the scalings associate with In terms of the 

exponential factor Bi, this dependence becomes exp ^— BjT ~g^ 3 ^ oc (T /10 8 K) S,//14 a Tf^ 1 g^ 1 ^ 28 ■ Since fjj does not 

scale strongly with To compared to the nuclear reaction rate (especially with heavier elements where Bi gets progressively 
larger), we ignore the scaling of fn with To and 514. To put this another way, the value of the number density, fjj, at the 
burning layer does not change as drastically as the value of the nuclear timescale at the burning layer with To. Putting 
these scalings all together, we have 

/ \ l+Si 

r coli =r coli , ^ 00gcm _ 2 j T e6 g u , (50) 

Table 1 lists the values of the scalings, the prefactor, calculated both numerically and analytically, and the largest effective 
temperature where DNB in diffusive equilibrium is valid. Above T 6i gDNB, equation (50) is not valid since DNB no longer 
occurs in diffusive equilibrium. 

5.3. Application to Cen X-4 in Quiescence 

The large number of transiently accreting NS also provide an astrophysical site for our work. The work of Rutledge et al. 
(1999) (sec Bildstcn and Rutledge 2001 for an overview) first showed that much of the quiescent emission from these objects 
was thermal emission from the surface, allowing for many T e measurements (Rutledge et al. 2001 a,b, 2000, 1999). Of all 
these LMXBs, only Cen X-4 is safely in the DNB regime. Its effective temperature is low enough (T e = 8.8 x 10 5 Rutledge 
et al. 2001a) so that DNB in diffusive equilibrium is a good approximation for a hydrogen column yn < 10 8 gcm~ 2 . The 
column depth of H on Cen X-4 after an outburst is expected to be yn < 10 8 gem -2 ; set by ignition conditions of type I 
X-ray bursts (Brown et al. 2002). 

Accretion outbursts were observed in 1969 and 1979 from Cen X-4, but the latter outburst was particularly weak. 
Presuming no outburst has occurred in the intervening years, this gives an age of 23 years. For a hydrogen column of 
jjh = 10 8 gem -2 after the 1979 outburst, we find that the present parameters of Cen X-4 would be yn ~ 2 x 10 6 gcm~ 2 
and T c « 4 x 10 7 K if there has been no accretion in quiescence. A much larger H column « 2 x 10 7 gcm~ 2 ) is needed 
to burn matter via DNB at the supplied rate if the NS is accreting in quiescence at a rate M q = 10 _13 M Q yr _1 . 

Presuming little accretion in quiescence and fixed T c , we evolve yn and the resulting flux as a function of the time since 
the last outburst. The age of the envelope or time since the last outburst, i ag0 is the age derived from the scaling of r co i 
with yn in the analytic solution, 

((dy H y l , 12 

ta.se = ) d VH = —T co \- (51) 

In Figure 8, we plot the resultant flux from a total integrated column of H, yn , against the lifetime of that column for 
fixed core temperature, T c . Evolution via DNB predicts that the outgoing flux will drop by ~ 12% over a twenty year 
period roughly 100 years after the outburst. This drop is due to the consumption of H that makes up the sensitivity strip 
of the H layer. At present (twenty years after the outburst) the flux will vary by only 3% over a ten year timescale. The 
associated Mdnb — 2x 10 _15 M Q yr _1 for this time after the outburst. Observations of Cen X-4 have placed constraints 
on the variation of T e over a five year timescale of < 10% (Rutledge et al. 2001a). Observational confirmation of DNB on 
Cen X-4 is therefore unlikely. 

6. CONCLUSIONS AND DISCUSSIONS 

We have shown that diffusive nuclear burning (DNB) is an effective mechanism that burns surface hydrogen on a NS 
on astrophysically relevant timescales. Our numerical and analytic solutions for the rates of DNB are in the limit where 
the H is in diffusive equilibrium with the underlying proton capturing elements. Our work is similar to that originally 
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carried out by Michaud and Fontaine (Michaud et al. 1984; Michaud & Fontaine 1984) for diffusion-induced hydrogen 
burning in white dwarfs and also by Iben & MacDonald (1985a, b). However, there are several crucial differences: the NS 
temperatures are much higher and the local gravity much stronger. Hence, we were able to make simplifying assumptions 
and avoid the active diffusion calculation with nuclear burning performed by Iben & MacDonald (1985a). We have also 
solved for the concentration profile in the degenerate regime self-consistently, whereas the effects of degeneracy were 
neglected in Michaud and Fontaine's paper. 

The largest remaining uncertainty in applying our work is the unknown composition of the material underlying the 
hydrogen. Our assumptions of low magnetic fields and surface temperatures (T e < 10 6 K for DNB to be nuclear-rate 
limited) limit the DNB applications to quiescent NSs in transient low-mass X-ray binaries and millisecond radio pulsars. 
Little is known about the surface composition of millisecond radio pulsars, though in the absence of our work one would 
definitely expect H to be dominant on the surface since these objects underwent mass accretion prior to becoming a 
pulsar. If there is enough carbon in the underlying composition from previous H/Hc burning during type I X-ray bursts 
(Schatz et al. 1999), then there is clearly adequate time to burn off the H. A clear detection of H at the photosphere of a 
millisecond radio pulsar would thus constrain the underlying abundances. 

The impact of an intervening helium layer depends on its thickness. A thin layer of helium, which does not penetrate 
into the burning layer and remains in the nondcgenerate regime, enhances the burning rate. This is because the electric 
field in a nondcgenerate helium plasma (E = 4m p g/3e) is smaller than the electric field in a nondcgenerate carbon plasma 
(E = I2m p g/7e). Hence the diffusive tail of hydrogen can penetrate more easily into the helium layer onto the underlying 
layer of proton capturing elements and the number fraction of hydrogen is enhanced at the burning layer. For a thick 
helium layer which goes beyond the burning layer, DNB would be effectively shut off. However, a thick helium layer would 
not likely sit on top of a thick proton capturing layer, if the proton capturing elements have the same A/Z as helium 
and both layers are degenerate. This is because the electric field is no longer a differentiating factor allowing the helium 
to mix downward with the proton capturing material. This would effectively dilute the abundance of proton capturing 
elements, thereby reducing the rate, but not shutting off DNB. The timescale for this mixing and its exact impact on the 
overall burning rate still needs more careful attention. 

The recent observations of a pair of X-ray spectral lines on 1E1207. 4-5209 (Sanwal et al. 2002; Mereghetti et al. 2002) 
show that hydrogen is not present on the surface of this young NS (Sanwal et al. 2002; Hailey & Mori 2002). These spectral 
lines appear to be mid-atomic elements like oxygen or neon (Hailey & Mori 2002) or helium if B ~ 10 14 G (Sanwal et al. 
2002). Since the age of the associated SN remnant is 7 kyrs, the mechanism of hydrogen removal must be extremely 
fast. However, our calculation is not directly applicable to this system. The spindown of 1E1207. 4-5209 implies a dipolc 
B-field of B « 3 x 10 12 G (Pavlov et al. 2002). This object also has a high temperature of T w 1.4 - 1.9 MK at which our 
assumption of diffusive equilibrium breaks down. However, the physics remains the same and we expect this mechanism 
to be active in this system. 

Since we have explored one special limit of DNB, it is not surprising that we have not found a physical system in 
which DNB is directly observable. In a future paper we will consider the problem of diffusion limited DNB, which is 
applicable to systems where the assumption of diffusive equilibrium breaks down (T e > 10 6 K for H/ 12 C). We will also 
consider the problem of DNB in highly magnetized sources such as young radio pulsars (B ~ 10 12 ~ 13 G) and magnetars 
(B ~ 10 14 - 15 G). 

We thank D. Uzdensky for showing us a nice mathematical trick, and R. Sunyacv for reminding us of the importance 
of the magnetic field on the thermal structure in radio pulsars. We would also like to thank S. W. Davis, C. J. Deloye, 
A. Socrates and D. Townsley for discussions and the referee for important clarifications. P.C. would like to thank the 
Department of Physics and Department of Astronomy at Columbia University, where the early and latter parts of this 
work was done, for their hospitality. This research was supported by NASA via grant NAG 5-8658 and by the NSF under 
Grants PHY99-07949 and AST01-96422. L. B. is a Cottrell Scholar of the Research Corporation. 
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H-layer 



H/C boundary 



C-layer \ 









A 



Diffusive H-tail into C-layer 

Fig. 1. — Diagram of a H/C envelope in diffusive equilibrium. The diffusive tail of hydrogen extends deep into the carbon, reaching 
temperatures where the hydrogen rapidly captures onto carbon. 



Table 1 

Power law exponents and prefactors and range of validity for equation (50). 



Reaction 


Igio ( r coi„o) [yrs] 




l + 6i 


T e fi DNB [K] 




numerical analytic 






12 C( P , 7 ) 13 N 


0.57 0.58 


136.93 


-5/12 


1 


14 N(p, 7 ) 15 


2.02 2.15 


152.31 


-6/14 


1.2 


16 0( P , 7 ) 17 F 


3.06 3.23 


166.96 


-7/16 


1.3 



12 



Chang, P. & Bildsten, L. 



p_ure_ C 




pure H 



1.5 2 2.5 3 

lg 10 y [g cm- 2 ] 

Fig. 2. — Electric field strength in an atmosphere in diffusive equilibrium. The electric field is shown for a total column of hydrogen of 
yil = 100 gem -2 and T e = 8 X 10 5 K. The electric field changes from the value for a pure hydrogen atmosphere to that for a pure carbon 
atmosphere. The bottom graph shows the mass fraction of the two components, Xf, as a function of column. The variation in Xi traces out 
the variation in the electric field. The electric field varies from one limit to another over a zone of order of a scale height, which agrees with 
the results of Alcock (1980). 
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Fig. 3. — Thermal structure and composition for a NS in diffusive equilibrium with T e = 8 X 10 5 K, 314 = 2, and yu = lOOgcm" 2 . 
The solid line is the numerical solution for the temperature and hydrogen number fraction, fjj = nfj/ntot, as a function of column depth. 
The dashed line is the approximate solution for the temperature and hydrogen number fraction, fjj. The approximate analytic solution for 
temperature is given by equation (16) and (24). The approximate solution for fu is given by equations (31) and (35) with = 1.5. 
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Fig. 4. — Differential hydrogen column burning rate taking into account p-p capture and p + 12 C capture. The bottom graph shows the 
number fraction (solid line) and temperature (dotted line). This model has yu = 100 gem -2 and T e = 8 X 10 5 K. The integrated burning 
rate for this model is yul T <io\ = 0.24 g cm -2 yr _1 . The burning peak occurs at a column of {/burn ~ 10 6 gcm -2 , where T = 2.9 X 10 7 K and 
p = 4.9 X 10 4 gem -3 . The local drift time ("n on = 10 4 s) is much shorter than the local burning time (th = 1-5 X 10 7 s). 
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Fig. 5. — Lifetime of a hydrogen column, r co i, and total mass burning rate, MdnBi as a function of the size of the hydrogen column for 
different fixed core temperatures and 314 = 2. For each model, we list the logarithmic core temperature and associated logarithmic effective 
temperature. The bottom plot shows the total mass burning rate for a NS with a fiducial radius of fOkm for a hydrogen column of yu- For 
a central temperature of lgioT c = 7.7, our assumption of constant flux breaks down for columns greater than 3 X f0 7 gcm~ 2 . At these large 
columns the heat release from nuclear burning becomes comparable to the flux. The power law dependence between the lifetime and yu is 
universal and the scaling does not change with different core temperatures. The derivation of this and other power law scalings is described 
in § 5.2. 
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Fig. 6. — Thermal structure of a NS with a total integrated H column, yu K, 100gcm~ 2 with different T e fi = (0.8, 1, 1.26) respectively. 
The burning layer for each of these models is represented by area between two vertical lines. Also plotted is the line Tj on = tjj, above which 
DNB docs not occur in diffusive equilibrium. 
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Fig. 7. — Comparison of scaling laws for the lifetime of the hydrogen column, r co i- The numerical solutions (solid lines) agree with a high 
degree of precision with the analytic scaling laws (dashed lines). The comparison for effective temperature and local gravity is given for a 
fixed hydrogen column of j/jj = 100 gem . For the plot of lifetime against total hydrogen column, the numerical calculation was done with 
a constant core temperature. The disagreement between analytic and numerical results for low effective temperatures is due to our chosen 
To,6 = 40 for the power law expansion of the burning rate. Expanding around a lower value of To,6 resolves this discrepancy. 
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